using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;

using Atomic.Interfaces.Atat;
using Atomic.Libraries.Jobs;
using Atomic.Libraries.Plotting;
using Atomic.Structures;
using Atomic.Thermodynamics.StateEquations;
using Atomic.Vasp;
using Atomic.Vasp.Options;

namespace Atomic.Samples.Structures
{
	public static class Structures1
	{
		public static void Sample1()
		{
			// Build a primitive hexagonal cell of iron atoms with nitrogen or vacancy at interstitial positions.
			Structure lattice0 = StructureHelper.CreateHexagonalStructureOctahedral(2.8, Atom.Fe, Atom.Placeholder(Atom.C, Atom.N, Atom.Vacancy));

			// Transform to conventional 6-atom cell by repeating the lattice vectors.
			Structure lattice = lattice0.Repeat(2, 1, 0, -1, 1, 0, 0, 0, 1);

			// Show the two structures.
			Structure.Show(lattice0, lattice);

			// Compute all symmetrically distinct combinations of filling the interstitial sites.
			Structure[] structures = lattice.FillPlaceholdersCombinatorialReduced().ToArray();

			// An example of a powerful Linq operation on the computed combinations.
			// First, vacancies are removed, and carbon and nitrogen atoms are counted and stored in a temporary anonymous type.
			// Then a filtering is performed.
			// Finally, sorting is done, first by carbon atoms, then by nitrogen atoms.
			Structure[] selectedStructures = structures
				.Select(s => new { Structure = s.RemoveVacancies(), C = s.Sites.Count(site => site.Atom == Atom.C), N = s.Sites.Count(site => site.Atom == Atom.N) })
				.Where(s => s.C + s.N <= 3)
				.OrderBy(s => s.C + s.N)
				.ThenBy(s => s.C)
				.Select(s => s.Structure)
				.ToArray();

			// Print number of structures and chemical formula for each of them.
			Console.WriteLine(selectedStructures.Length);
			foreach (Structure structure in selectedStructures)
			{
				Console.WriteLine(structure.Formula);
			}

			// Visualize the selected structures.
			Structure.Show(selectedStructures);
		}

		public static void Sample2()
		{
			// Create a fcc structure with a fixed Fe and a site with C/N/Vac.
			Structure s0 = Structure.CreateFaceCenteredCubicStructure(4.0)
				.Add(Atom.Fe, 0.0, 0.0, 0.0)
				.Add(Atom.Placeholder(Atom.C, Atom.N, Atom.Vacancy), 0.5, 0.5, 0.5)
				.Repeat(2, 1, 1);

			// Iterate through all symmetrically distinct combinations.
			foreach (Structure s in s0.FillPlaceholdersCombinatorialReduced())
			{
				// ...
			}

			// Or just visualize them all. Use LINQ here to strip all vacancies.
			Structure.Show(s0.FillPlaceholdersCombinatorialReduced().Select(s => s.RemoveVacancies()));
		}

		public static void Sample3()
		{
			// Create fcc structure with two atoms per cell.
			Structure s = Structure.CreateFaceCenteredCubicStructure(4.0)
				.Add(Atom.Fe.SpinPolarize(3.0), 0.0, 0.0, 0.0)
				.Add(Atom.C.SpinPolarize(0.5), 0.5, 0.5, 0.5);

			// Prepare VASP options represented as options.
			VaspProblem vasp = new VaspProblem();
			vasp.Structure = s;
			vasp.Convergence.PlaneWaveCutOff = 350.0;
			vasp.Convergence.EnergyBreakCondition = 1.0e-5;
			vasp.Convergence.Precision = VaspPrecision.Accurate;
			vasp.Options.Add(new VolumeIonicRelaxation(40));
			vasp.Options.Add(new MethfesselPaxtonSmearing(1, 0.2));
			vasp.Options.Add(new MonkhorstPackMesh(11, 11, 11));

			// Start VASP and extract the result.
			VaspResult result = vasp.Run();

			// Extract volume and lattice constant of the optimized structure.
			double v = result.FinalStructure.Volume;
			double a =
				((FaceCenteredCubicLattice)result.FinalStructure.Lattice).A0;

			Console.WriteLine(v);
			Console.WriteLine(a);
		}
	
		public static void Sample4()
		{
			// Create a base structure with one atom.
			Structure s0 = Structure.CreateFaceCenteredCubicStructure(4.04)
				.Add(Atom.Al, 0.0, 0.0, 0.0);

			List<VolumeEnergyPoint> points = new List<VolumeEnergyPoint>();

			foreach (double a in new double[] { 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04 })
			{
				// Scale structure linearly.
				Structure s = s0.Scale(a);

				// Prepare VASP options represented as objects.
				VaspProblem vasp = new VaspProblem();
				vasp.Structure = s;
				vasp.Options.Add(new MethfesselPaxtonSmearing());
				vasp.Options.Add(new MonkhorstPackMesh(11, 11, 11));

				// Run VASP. An object is generated containing relevant output.
				VaspResult result = vasp.Run();

				points.Add(new VolumeEnergyPoint(result.FinalStructure.Volume, result.Energy));
			}

			// Fit an equation of state.
			IStateEquation bm4 = BirchMurnaghan4.Fit(points);

			// Plot with the original data points.
			Gnuplot.Plot("Volume, Å^3", "Energy, eV",
				new DataPlot(points),
				new FunctionPlot(x => bm4.Energy(x), 15.0, 18.6, new PlotLabel("BM4")));
		}

		public static void Sample5_Part1(string[] args)
		{
			// This is the entry point of submitted jobs. Keep this as the first part of the program.
			if (JobManager.TryParseCommandLine(args))
			{
				// Execution was handled by the job manager. Just exit now.
				return;
			}

			// Create a base structure with one atom.
			Structure s0 = Structure.CreateFaceCenteredCubicStructure(4.04)
				.Add(Atom.Al, 0.0, 0.0, 0.0);

			List<Job> jobs = new List<Job>();

			foreach (double a in new double[] { 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04 })
			{
				// Scale structure linearly.
				Structure s = s0.Scale(a);

				// Prepare VASP options represented as objects.
				VaspProblem vasp = new VaspProblem();
				vasp.Structure = s;
				vasp.Options.Add(new MethfesselPaxtonSmearing());
				vasp.Options.Add(new MonkhorstPackMesh(11, 11, 11));

				// Don't start VASP, just add the job to the list. Request 1 hour on 1 processor.
				jobs.Add(new VaspJob(vasp, a.ToString(), new JobResources(TimeSpan.FromHours(1.0), 1)));
			}

			// Submit all jobs.
			JobManager.Submit(jobs, "atomic");
		}

		public static void Sample5_Part2()
		{
			// Load previously computed jobs.
			JobManager manager = JobManager.Load(new DirectoryInfo("/home/mbjba/cluster/out"));

			List<VolumeEnergyPoint> points = new List<VolumeEnergyPoint>();

			foreach (VaspJob job in manager.Jobs)
			{
				// Retrieve output from VASP.
				VaspResult result = job.Result;

				points.Add(new VolumeEnergyPoint(result.FinalStructure.Volume, result.Energy));
			}

			/*// Fit an equation of state.
			BirchMurnaghan5 bm5 = BirchMurnaghan5.Fit(points);

			// Plot with the original data points.
			Gnuplot.Plot(new AxisLabel("Volume, Å^3"), new AxisLabel("Energy, eV"),
				new DataPlot(points),
				new FunctionPlot(bm5, 15.0, 18.6, new PlotTitle("BM5")));*/
		}
	}
}
